Tesis

#Instalación rSFFreader
# source("http://bioconductor.org/biocLite.R")
# biocLite("rSFFreader")
require(rSFFreader)
library(gdata)
library(reshape2)
library(plotly)
library(ggplot2)
library(Biostrings)
#source("https://bioconductor.org/biocLite.R")
#biocLite("R453Plus1Toolbox")
library(R453Plus1Toolbox)
samples<-grep(".sff",list.files("./samples"),value=TRUE)

Ejemplo de libro

#ejemplo de libro
sff<- load454SampleData()

Total number of reads to be read: 1000
reading header for sff file:/home/sergio/R/x86_64-pc-linux-gnu-library/3.3/rSFFreader/extdata/Small454Test.sff
reading file:/home/sergio/R/x86_64-pc-linux-gnu-library/3.3/rSFFreader/extdata/Small454Test.sff

##Generate some QA plots:
##Read length histograms:
par(mfrow=c(2,2))
clipMode(sff) <- "raw"
hist(width(sff),breaks=500,col="grey",xlab="Read Length",main="Raw Read Length")
## Base by position plots:
clipMode(sff) <- "raw"
ac <- alphabetByCycle(sread(sff),alphabet=c("A","C","T","G","N"))
ac.reads <- apply(ac,2,sum)
acf <- sweep(ac,MARGIN=2,FUN="/",STATS=apply(ac,2,sum))
matplot(cbind(t(acf),ac.reads/ac.reads[1]),col=c("green","blue","black","red","darkgrey","purple"),
        type="l",lty=1,xlab="Base Position",ylab="Base Frequency",
        main="Base by position")

cols <- c("green","blue","black","red","darkgrey","purple")
leg <- c("A","C","T","G","N","%reads")
legend("topright", col=cols, legend=leg, pch=18, cex=.8)
clipMode(sff) <- "full"
hist(width(sff),breaks=500,col="grey",xlab="Read Length",main="Clipped Read Length")

ac <- alphabetByCycle(sread(sff),alphabet=c("A","C","T","G","N"))
ac.reads <- apply(ac,2,sum)
acf <- sweep(ac,MARGIN=2,FUN="/",STATS=apply(ac,2,sum))
matplot(cbind(t(acf),ac.reads/ac.reads[1]),col=c("green","blue","black","red","darkgrey","purple"),
        type="l",lty=1,xlab="Base Position",ylab="Base Frequency",
        main="Base by position")
legend("topright", col=cols, legend=leg, pch=18, cex=.8)

IMPORT DATA

SIN CLIP - RAW DATA ANALYSIS

lenStat_raw<-vector()
lenStat_qual<-vector()
nucFreq_raw<-matrix(nrow=1000, ncol=1)
nucFreq_qual<-matrix(nrow=1000, ncol=1)
count_raw<-data.frame()
count_qual<-data.frame()
len_raw<-list()
len_qual<-list()
for(i in samples){
  sff<-readSff(paste("samples/",i,sep=""), use.qualities=TRUE, use.names=TRUE,clipMode = c("raw"), verbose=TRUE)
##Generate some QA plots:
  ##Read length histograms (with and without clipping):
  
  # RAW
  par(mfrow=c(2,2))
  clipMode(sff) <- "raw"
  hist(width(sff),breaks=500,col="grey",xlab="Read Length",
       xlim= c(50,800), main= "RAW read length")
  lenStat_raw<-rbind(lenStat_raw, c(gsub("454Reads.MID_","",gsub(".sff","",i)), summary(width(sff)) ) )
  count_raw[i,1]<-length(sff)
  len_raw[[i]]<-width(sff)
    ## Base by position plots:
  ac <- alphabetByCycle(sread(sff),alphabet=c("A","C","T","G","N"))
  ac.reads <- apply(ac,2,sum)
  acf <- sweep(ac,MARGIN=2,FUN="/",STATS=apply(ac,2,sum))
  tacf<-t(acf); colnames(tacf)<-paste(gsub("454Reads.MID_","",gsub(".sff","",i)),"_",  c("A","C","T","G","N"),sep="")
  nucFreq_raw<-cbindX(nucFreq_raw,tacf)
  matplot(cbind(t(acf),ac.reads/ac.reads[1]),col=c("green","blue","black","red","darkgrey","purple"),
        type="l",lty=1,xlab="Base Position",ylab="Base Frequency",
        main="Base by position", xlim=c(0,1500))
  cols <- c("green","blue","black","red","darkgrey","purple")
  leg <- c("A","C","T","G","N","%reads")
  legend("topright", col=cols, legend=leg, pch=18, cex=.8)
  
  ### FILTER QUALITY
  clipMode(sff) <- "full"
  hist(width(sff),breaks=500,col="grey",xlab="Read Length",
       xlim= c(50,800),main="CLIPPED read length" )
  lenStat_qual<-rbind(lenStat_qual,c(gsub("454Reads.MID_","",gsub(".sff","",i)), summary(width(sff)) ) )
  count_qual[i,1]<-length(sff)
  len_qual[[i]]<-width(sff)
  ## Base by position plots:
  ac <- alphabetByCycle(sread(sff),alphabet=c("A","C","T","G","N"))
  ac.reads <- apply(ac,2,sum)
  acf <- sweep(ac,MARGIN=2,FUN="/",STATS=apply(ac,2,sum))
    tacf<-t(acf); colnames(tacf)<-paste(gsub("454Reads.MID_","",gsub(".sff","",i)),"_",  c("A","C","T","G","N"),sep="")
  nucFreq_qual<-cbindX(nucFreq_qual,tacf)
  matplot(cbind(t(acf),ac.reads/ac.reads[1]),col=c("green","blue","black","red","darkgrey","purple"),
        type="l",lty=1,xlab="Base Position",ylab="Base Frequency",
        main="Base by position", xlim=c(0,1000))
  par(mfrow=c(1,1))
  legend("topright", col=cols, legend=leg, pch=18, cex=.8)
  title(paste("sample: ",gsub("454Reads.MID_","",gsub(".sff","",i)),sep=""))
  dev.copy(png,filename=paste("Explo_",gsub("454Reads.MID_","",gsub(".sff","",i)),".png",sep=""));
  dev.off ();
}

Total number of reads to be read: 2233
reading header for sff file:samples/454Reads.MID_s01_rs3.sff
reading file:samples/454Reads.MID_s01_rs3.sff
Total number of reads to be read: 2405
reading header for sff file:samples/454Reads.MID_s02_rs3.sff
reading file:samples/454Reads.MID_s02_rs3.sff

Total number of reads to be read: 2054
reading header for sff file:samples/454Reads.MID_s03_rs5.sff
reading file:samples/454Reads.MID_s03_rs5.sff

Total number of reads to be read: 11821
reading header for sff file:samples/454Reads.MID_s04_rl3.sff
reading file:samples/454Reads.MID_s04_rl3.sff

Total number of reads to be read: 10580
reading header for sff file:samples/454Reads.MID_s05_rl3.sff
reading file:samples/454Reads.MID_s05_rl3.sff

Total number of reads to be read: 8869
reading header for sff file:samples/454Reads.MID_s06_rl5.sff
reading file:samples/454Reads.MID_s06_rl5.sff

Total number of reads to be read: 9255
reading header for sff file:samples/454Reads.MID_s07_dl3.sff
reading file:samples/454Reads.MID_s07_dl3.sff

Total number of reads to be read: 9579
reading header for sff file:samples/454Reads.MID_s08_dl3.sff
reading file:samples/454Reads.MID_s08_dl3.sff

Total number of reads to be read: 9987
reading header for sff file:samples/454Reads.MID_s09_dl5.sff
reading file:samples/454Reads.MID_s09_dl5.sff

Estadísticas de las lecturas

#Statistics about lenght
as.data.frame(lenStat_raw)

write.csv(as.data.frame(lenStat_raw), "Raw_Readlength_stats.csv",row.names = FALSE )
as.data.frame(lenStat_qual)

write.csv(as.data.frame(lenStat_qual), "Clip_Readlength_stats.csv",row.names = FALSE )
#Nucleotide frequency by position
write.csv(nucFreq_raw, "Raw_NuclFreqByPos.csv",row.names = FALSE )
write.csv(nucFreq_qual, "Clip_NuclFreqByPos.csv",row.names = FALSE )
#All length data by sample
count_raw$sample<-gsub("454Reads.MID_","",gsub(".sff","",rownames(count_raw)))
count_qual$sample<-gsub("454Reads.MID_","",gsub(".sff","",rownames(count_qual)))
length_raw<-do.call(cbindX, lapply(len_raw, as.data.frame))
colnames(length_raw)<-gsub("454Reads.MID_","",gsub(".sff","",rownames(count_raw)))
length_qual<-do.call(cbindX, lapply(len_qual, as.data.frame))
colnames(length_qual)<-gsub("454Reads.MID_","",gsub(".sff","",rownames(count_qual)))
rownames(count_raw)<-NULL;rownames(count_qual)<-NULL
write.csv(length_raw, "Raw_AllReadlength.csv",row.names = FALSE )
write.csv(length_qual, "Clip_AllReadlength.csv",row.names = FALSE )
#AMount of sequences by sample
count_raw

write.csv(count_raw, "Raw_CountBySample.csv",row.names = FALSE )
count_qual

write.csv(count_qual, "Clip_CountBySample.csv",row.names = FALSE )

#GRAFICO DE BARRAS CON LA CANTIDAD DE LECTURAS
#plot(cantidad$cantidad, type="b")
library(ggplot2)
p<-ggplot(count_raw, aes(x=sample, weight=V1))+ geom_bar(fill="#2b8cbe")+
  theme(axis.text.x = element_text(angle = 45, hjust = 1))
ggplotly(p)

#GRAFICO DE CAJAS
#boxplot(longitud)
p<-ggplot(data = melt(length_raw), aes(x=variable, y=value)) + geom_boxplot(aes(fill=variable))+
  theme(axis.text.x = element_text(angle = 45, hjust = 1))+ ggtitle("Raw Reads")
ggplotly(p)

#GRAFICO DE CAJAS
#boxplot(longitud)
p<-ggplot(data = melt(length_qual), aes(x=variable, y=value)) + geom_boxplot(aes(fill=variable))+
  theme(axis.text.x = element_text(angle = 45, hjust = 1))+ ggtitle("Clipped Reads")
ggplotly(p)

for(i in samples){
  sff<-readSff(paste("samples/",i,sep=""), use.qualities=TRUE, use.names=TRUE,clipMode = c("raw"), verbose=TRUE)
customClip(sff) <- IRanges(start = 1, end = 15)
clipMode(sff) <- "custom"
print("distinct MIDs at 5'")
print(length(table(counts=as.character(sread(sff)))))
print("More frequents")
print(sort(table(counts=as.character(sread(sff))), decreasing=TRUE))
}

Total number of reads to be read: 2233
reading header for sff file:samples/454Reads.MID_s01_rs3.sff
reading file:samples/454Reads.MID_s01_rs3.sff
[1] "distinct MIDs at 5'"
[1] 3
[1] "More frequents"
counts
TCAGTGTACTACTCC TCAGTGTACTACTCT TCAGTGTACTACTCA 
           2224               5               4 
Total number of reads to be read: 2405
reading header for sff file:samples/454Reads.MID_s02_rs3.sff
reading file:samples/454Reads.MID_s02_rs3.sff
[1] "distinct MIDs at 5'"
[1] 4
[1] "More frequents"
counts
TCAGACGACTACAGC TCAGACGACTACAGT TCAGACGACTACAGG TCAGACGACTACAGA 
           2396               4               3               2 
Total number of reads to be read: 2054
reading header for sff file:samples/454Reads.MID_s03_rs5.sff
reading file:samples/454Reads.MID_s03_rs5.sff
[1] "distinct MIDs at 5'"
[1] 2
[1] "More frequents"
counts
TCAGAGCACTGTAGC TCAGAGCACTGTAGT 
           2053               1 
Total number of reads to be read: 11821
reading header for sff file:samples/454Reads.MID_s04_rl3.sff
reading file:samples/454Reads.MID_s04_rl3.sff
[1] "distinct MIDs at 5'"
[1] 3
[1] "More frequents"
counts
TCAGCACACGATAGC TCAGCACACGATAGA TCAGCACACGATAGT 
          11794              20               7 
Total number of reads to be read: 10580
reading header for sff file:samples/454Reads.MID_s05_rl3.sff
reading file:samples/454Reads.MID_s05_rl3.sff
[1] "distinct MIDs at 5'"
[1] 4
[1] "More frequents"
counts
TCAGCACTCGCACGC TCAGCACTCGCACGA TCAGCACTCGCACGT TCAGCACTCGCACGG 
          10559              12               7               2 
Total number of reads to be read: 8869
reading header for sff file:samples/454Reads.MID_s06_rl5.sff
reading file:samples/454Reads.MID_s06_rl5.sff
[1] "distinct MIDs at 5'"
[1] 3
[1] "More frequents"
counts
TCAGCAGACGTCTGC TCAGCAGACGTCTGA TCAGCAGACGTCTGT 
           8853              12               4 
Total number of reads to be read: 9255
reading header for sff file:samples/454Reads.MID_s07_dl3.sff
reading file:samples/454Reads.MID_s07_dl3.sff
[1] "distinct MIDs at 5'"
[1] 3
[1] "More frequents"
counts
TCAGATCGTCTGTGC TCAGATCGTCTGTGA TCAGATCGTCTGTGT 
           9228              18               9 
Total number of reads to be read: 9579
reading header for sff file:samples/454Reads.MID_s08_dl3.sff
reading file:samples/454Reads.MID_s08_dl3.sff
[1] "distinct MIDs at 5'"
[1] 4
[1] "More frequents"
counts
TCAGATGTACGATGC TCAGATGTACGATGA TCAGATGTACGATGT TCAGATGTACGATGG 
           9562              11               5               1 
Total number of reads to be read: 9987
reading header for sff file:samples/454Reads.MID_s09_dl5.sff
reading file:samples/454Reads.MID_s09_dl5.sff
[1] "distinct MIDs at 5'"
[1] 4
[1] "More frequents"
counts
TCAGATGTGTCTAGC TCAGATGTGTCTAGA TCAGATGTGTCTAGT TCAGATGTGTCTAGG 
           9969              13               4               1 

Explore Each seque

quantile(as(readFullQ, "numeric"))

  0%  25%  50%  75% 100% 
   0   11   13   25   40 

readsRaw<-sread(sff)
readsRawQ<-   as(quality(sff), "PhredQuality")  
library(ShortRead); setSR <- ShortReadQ(sread=readsRaw, quality=FastqQuality(BStringSet(readsRawQ)), BStringSet(readsRawQ))
myMA <- as(quality(setSR), "matrix")

Processed with SSFinfo

file_name<-gsub("454Reads.MID_","",gsub(".sff","",i))
library(Biostrings)
pro1_read = readDNAStringSet(paste("samples/processed_files/",file_name,"/",file_name,".fna", sep=""), format="fasta")
pro1_read_width<-as.data.frame(cbind(names(pro1_read),width(pro1_read)))
head(pro1_read_width)

sffContainer <- readSFF(paste("samples/",i,sep=""))

Reading file 454Reads.MID_s01_rs3.sff ... done! 

showClass("SFFContainer")

Class "SFFContainer" [package "R453Plus1Toolbox"]

Slots:
                                                                                    
Name:                       name            flowgramFormat                 flowChars
Class:                 character                   numeric                 character
                                                                                    
Name:                keySequence           clipQualityLeft          clipQualityRight
Class:                 character                   numeric                   numeric
                                                                                    
Name:            clipAdapterLeft          clipAdapterRight                 flowgrams
Class:                   numeric                   numeric                      list
                                                          
Name:                flowIndexes                     reads
Class:                      list QualityScaledDNAStringSet

reads(sffContainer)

  A QualityScaledDNAStringSet instance containing:

  A DNAStringSet instance of length 2233
       width seq                                                                      names               
   [1]   684 TCAGTGTACTACTCCACGACGTTTGTAAACGACGC...TACGTTAGGGTAAACGGAACCGGACGGAAACGGG H3C8HXX01ATG0S
   [2]   552 TCAGTGTACTACTCCACGACGTTTGTAAAACGACG...AACAAGCATATTGACGCTGGGAAAGGACCACGAG H3C8HXX01BDS7I
   [3]   546 TCAGTGTACTACTCCACGACGTTTGTAAAACGACG...CAAAAACAGCATATGACGCTGGGAAGACCGAGAG H3C8HXX01ANWA5
   [4]   769 TCAGTGTACTACTCCACGACGTTTGTAAACGACGC...AAAACGGGAAACCCGGTAACGGGAAGGGAAACCG H3C8HXX01A3BQX
   [5]   564 TCAGTGTACTACTCCACGACGTTTGTAAAACGACG...AAGAACATATATACGTCNCGGAAGAACCGAAGAG H3C8HXX01AMUIH
   ...   ... ...
[2229]   580 TCAGTGTACTACTCCACGACGTTTGTAAAACGACG...ANTTTAGCGTACTGGGTAAAAGGACCCTAGACGG H3C8HXX01AQGHL
[2230]   700 TCAGTGTACTACTCCACGACGTTTGTAAAACGACG...AGGACCTACGGGTAAAGGAAACCGTACCGGGACG H3C8HXX01BBPGN
[2231]   616 TCAGTGTACTACTCCACGACGTTTGTAAAACGACG...GGTAACCGNGTTGGGAAAACGGACCGTAGAACGG H3C8HXX01AR5XF
[2232]   802 TCAGTGTACTACTCCACGACGTTTGTAAAACGACG...CGACCGGGTTAAAAGGACCCGGGTAAAGTACCGG H3C8HXX01BBAJR
[2233]   566 TCAGTGTACTACTCCACGACGTTGTAAAACGAACG...AATTGTACGATTCCGGGGAAAGGAAACCGACGAG H3C8HXX01A1VEZ

  A PhredQuality instance of length 2233
       width seq                                                                      names               
   [1]   684 FFFFFFFFFFFFFFFIIIIIIHHHIIHHHIIIIII...,.0..,,,,,,,,.0,,,,,,,,,,,,...,,,, H3C8HXX01ATG0S
   [2]   552 FFFFFFFFFFFFFFFIIIIHHDDDFE::::EFFHH...,,,)),003111444..,,,,,,,,,,,,,,0.. H3C8HXX01BDS7I
   [3]   546 IIIIIIIIIIIIIIIIIIIIIHHHIICCCAIIFII...(,-//.4+1121//115/,,,,,,,,.,,,,,// H3C8HXX01ANWA5
   [4]   769 FFFFFFFFFFFFFFFIIIIIIHHHIIHHHIIIIIH...,,,,.*))))).,,,,,,,,,,,,,,,-//2770 H3C8HXX01A3BQX
   [5]   564 FFFFFFFFFFFFFFFIIIIIIHFHHH;:::FFFFF...))6002.,,,,13,,,,!,,,,,,,,,,....,, H3C8HXX01AMUIH
   ...   ... ...
[2229]   580 FFFFCCCFFDDDDDDEFFHHIHD@BB@@AAFD<::...,!,,,,,,,,,,,0,,,,,,,,,),,,,.0,,,, H3C8HXX01AQGHL
[2230]   700 FFFFFFFFD@@@FFFHHFFFGFCCCA;;;;::4=A...,,,,,00,,,,0./..,,,,030,,,,,,,,,,, H3C8HXX01BBPGN
[2231]   616 DD?:::ADDFFFFFFFFFFGFF??DA;;999?6?<...00.,,,,.!,,,,,/..,,,,,),,20,,,,,,, H3C8HXX01AR5XF
[2232]   802 FFFFFFFFFFDDFAAAAA:::;88822222//-//....3444/,,,,,,,,,,),,,,,,,,,.-4-,,,, H3C8HXX01BBAJR
[2233]   566 DFFCCCC???FDDAA?AAB=555544464?=55=?...,,,,,,04.---882222444---,,,,,,,,,, H3C8HXX01A1VEZ

#qualityReportSFF(sffContainer, "report.pdf")
positionQualityBoxplot(sffContainer)

#dinucleotideOddsRatio
#sff2fastq()

#Example flowgram
seqtab<-cbindX(as.data.frame(sffContainer@flowgrams$H3C8HXX01ATG0S),
as.data.frame(sffContainer@flowIndexes$H3C8HXX01ATG0S),
as.data.frame(sffContainer@reads$H3C8HXX01ATG0S),
as.data.frame(sffContainer@reads@quality$H3C8HXX01ATG0S),
as.data.frame(as.numeric(sffContainer@reads@quality$H3C8HXX01ATG0S)) )
colnames(seqtab)<-c("flowgrams","flowIndexes","basecall","quality_char","quality_value")
seqtab%>%
  group_by(basecall)%>%
        plot_ly( y = ~quality_value, color = ~basecall, type = "box")

seqtab%>%
  group_by(basecall)%>%
        plot_ly( y = ~flowgrams, color = ~basecall, type = "box")
  
     

#sffContainer@clipQualityLeft
#sffContainer@clipQualityRight


###
#convert http://www.drive5.com/usearch/manual/quality_score.html
#http://gatkforums.broadinstitute.org/gatk/discussion/4260/phred-scaled-quality-scores

#as.numeric(sffContainer@reads@quality$H3C8HXX01ATG0S)

#### Hacer una  vector para cada secuencia con todas las frecuencias de 4 nucleótidos

reads = sread(sff)
nf = oligonucleotideFrequency(reads, width=4)
#hclust(dist(nf)) # do hierarchical clustering of your tetra freq.